Emplacement of the Franklin large igneous province and initiation of the Sturtian Snowball Earth

During the Cryogenian (720 to 635 Ma ago) Snowball Earth glaciations, ice extended to sea level near the equator. The cause of this catastrophic failure of Earth’s thermostat has been unclear, but previous geochronology has suggested a rough coincidence of glacial onset with one of the largest magmatic episodes in the geological record, the Franklin large igneous province. U-Pb geochronology on zircon and baddeleyite from sills associated with the paleo-equatorial Franklin large igneous province in Arctic Canada record rapid emplacement between 719.86 ± 0.21 and 718.61 ± 0.30 Ma ago, 0.9 to 1.6 Ma before the onset of widespread glaciation. Geologic observations and (U-Th)/He dates on Franklin sills are compatible with major post–Franklin exhumation, possibly due to development of mafic volcanic highlands on windward equatorial Laurentia and increased global weatherability. After a transient magmatic CO2 flux, long-term carbon sequestration associated with increased weatherability could have nudged Earth over the threshold for runaway ice-albedo feedback.

20-30 min on a hot plate and then sonicated for the same amount of time before being rinsed and loaded into microcapsules for dissolution with ET535 spike. U-Pb dates and uncertainties were calculated following (79). The value 137.818 ± 0.045 (2σ) was used for the 238 U/ 235 U ratio in natural zircon (75). Model Th/U ratios were calculated iteratively from measured 206 Pb/ 208 Pb ratios and calculated 206 Pb/ 238 U ages. Up to 1 pg of common Pb (Pbc) was assumed to be procedural blank and accounted for using the measured laboratory Pbc isotopic composition. Excess Pbc was attributed to initial common Pb using the two-stage Pb isotope evolution model (80) at the nominal sample age. Concordia and upper-intercept dates were plotted using IsoplotR (71).
Sm-Nd isotope geochemistry methods Sm and Nd isotope geochemistry was done at the Isotope Geology Laboratory at Boise State University. A hundred mg of sample powder were spiked with a mixed 149 Sm-150 Nd tracer, dissolved with 1.8 mL 29M HF + 0.2 mL 15M HNO3 in Parr pressure vessels at 220°C for 18 hours, dried and redissolved in 0.5 mL concentrated HNO3 twice, then dried and redissolved in 2 mL 6M HCl at 180°C for 12 hours. Total dissolutions were dried and redissolved in 5 mL 1M HCl + 0.1M HF at 180°C overnight. Bulk rare earth elements were separated by standard dilute HCl and HNO3 based cation exchange chemistry on 6 mm internal diameter (i.d.) x 20 cm columns of AG-50W-X8 resin, H+ form, 200-400 mesh); Sm and Nd were separated by reverse phase HDEHP chromatography on 4 mm i.d. x 10 cm columns of Eichrom Ln-spec resin, 50-100 mesh, using the methods developed by (81). Sm and Nd isotopes were measured on a IsotopX Phoenix X62 TIMS in static and dynamic Faraday modes, respectively. Instrumental mass fractionation of Sm and Nd isotopes was corrected with an exponential law relative to 146 Nd/ 144 Nd = 0.7219 and 152 Sm/ 147 Sm = 1.7831. 143 Nd/ 144 Nd ratio is reported as spike-stripped and bias-corrected relative to the accepted value of JNdi-1 standard (0.512106).
(U-Th)/He Thermochronology (U-Th)/He thermochronology uses the radioactive decay of uranium (U), thorium (Th) and, to a much lesser extent, samarium (Sm) to helium (He) and the dependence of He diffusion behavior in minerals on temperature. A (U-Th)/He date represents the time-integrated history of He production and He loss over the mineral's time-temperature (t-T) history (e.g., 82). The thermal history can be related to geological processes such as tectonic exhumation, erosion, or burial when interpreted in the context of other geologic observations such as unconformable contacts and petrological temperature indicators. In addition to temperature, He diffusion in natural materials is dependent on physical characteristics of the mineral, such as parent-nuclide zonation, grain size, and degree of crystal lattice damage from recoil during radioactive decay ("radiation damage"). Of these, radiation damage has the largest systematic effect and for some thermal histories can result in different (U-Th)/He dates among grains from the same sample. This radiation damage dependence can be used to gain information about multiple parts of the sample's thermal history that correspond to different temperatures recorded by the different sensitivities of more or less damaged grains (e.g., 83). Radiation damage is quantified using the proxy "effective uranium" (eU) based on parent nuclide concentration as: eU = [U] + 0.238[Th] + 0.0012 [Sm] Characteristic trends in date-eU plots can aid interpretation of datasets and can be simulated using radiation damage accumulation and annealing models for apatite (e.g., RDAAM, 83) and zircon (e.g., ZRDAAM, 84).
Data Collection (U-Th)/He analysis was done at the University of Colorado Boulder Thermochronology Research and Instrumentation Lab (CU TRaIL). Individual zircon and apatite grains from sample F1966 were analyzed. The number of grains available for analysis was limited by the small size of the original sample. Zircon from other samples were unavailable for thermochronology due to annealing during geochronology data collection. AHe dates for additional samples were not obtained because apatite is more easily thermally reset than zircon, which for this region results in maximum AHe dates in the Phanerozoic (85)-inapplicable to the Neoproterozoic thermal history. Grains were selected using a Leica M614 binocular microscope, measured, characterized, and then packed in niobium tubes following standard procedures. Packed grains were then degassed for He measurement in an ASI Alphachron 774 extraction and measurement line. Degassed grains were dissolved and U, Th, and Sm measurements from dissolved solutions were made with an Agilent 7900 quadrupole inductively-coupled plasma mass spectrometer. Additional details of grain selection and analytical procedures for apatite and zircon are described in (86) and (87), respectively. Uncorrected dates were calculated from measured U, Th, Sm and He values and corrected for alpha-ejection based on mineral geometry following (88). Uncertainty on eU is conservatively assumed to be 15% of the calculated value. Date uncertainty is reported as the 2σ propagated analytical uncertainty of the U, Th, Sm, and He measurements. All data are reported in Table S6.

Inverse Thermal History Modeling and Hypothesis Testing
We carried out inverse thermal history modeling to test the consistency of the (U-Th)/He data with substantial exhumation and chemical weathering of the Franklin LIP shortly after its emplacement due to development of a mafic volcanic highland. (U-Th)/He dates were modeled using the HeFTy software (89) implementing the RDAAM and the ZRDAAM. Dates were grouped together by eU and averaged, which is a commonly used approach (90; Table S7). The two highest eU zircon grains were excluded from modeling due to greater uncertainty in He diffusion behavior at higher damage levels (84,91,92). Modeling returns possible t-T paths that fit the input data within prescribed goodness of fit criteria: >0.5 for all data to be considered a "good fit" and >0.95 for all data to be considered an "acceptable fit" (93). We used a hypothesistesting strategy for the two t-T simulations presented here and include key geologic constraints on the sample histories in the model frameworks (e.g., 90).
Both models start at the time of emplacement and require surface temperatures at 515 Ma in keeping with the regional Cambrian unconformity constraint. Phanerozoic constraints on the thermal history are drawn from existing literature, particularly (85,94,95), and are derived in part from estimates of former sedimentary cover from kimberlite sedimentary xenolith suites and other geologic observations (Table S7). The Phanerozoic portion of the path is tightly constrained by the AHe data reported here, particularly from 145 Ma to the present. The differences between models are therefore in the period between emplacement and 515 Ma.
Model A was designed to test compatibility of the observed (U-Th)/He dates with the hypothesis of Franklin sill emplacement at <5 km depth followed by rapid cooling related to exhumation shortly after emplacement-the "Neoproterozoic Exhumation Hypothesis" (Fig. S5A). The starting depth of <5 km was largely based on the geometry of the Mackenzie-age Fortress Gabbro sills, which the Douglas Peninsula sill was formerly considered a part of. Although the geometry of the Douglas Peninsula sill is not fully known, the surrounding Fortress Gabbro sills have been mapped as "nested cone-sheets" based on their saucer-shape morphologies (96), which requires emplacement at shallow depths of 1-5 km (97). The model starts at 110-60°C to reflect emplacement at a typical depth (2-4 km, temperature calculated assuming an ~25°C/km geotherm and an average surface temperature of 10°C), with no constraints imposed between emplacement and surface conditions at 515 Ma. This model yields numerous good-fit paths, including those that display major cooling and erosion shortly after emplacement. This outcome indicates that the data are consistent with the hypothesis of multiple kilometers of weathering and erosion shortly after sample emplacement, possibly due, as we argue, to development and weathering of mafic volcanic highlands.
An alternative model, Model B, the "Neoproterozoic Burial/Exhumation Hypothesis," allows for a broader range of emplacement temperatures and corresponding depths from 0-6 km followed by Neoproterozoic burial prior to exhumation by 515 Ma (rather than only Neoproterozoic exhumation). We ran this model to test the uniqueness of the t-T paths generated in Model A. Similar to Model A, this simulation produces many good-fit t-T paths (Fig. S5B). This outcome indicates that the (U-Th)/He data alone do not require substantial exhumation immediately after Franklin emplacement. This can also be seen in forward model predictions using representative paths from each model simulation (Fig. S5C). However, based on the larger geologic context as described in this paper, we consider the exhumation scenario of Model A to be more geologically likely.   (27) for comparison. Most samples in this study plot similarly to Type 2 but S8 and 17RAT-R35B1 show evidence for significant crustal contamination and assimilation, more similar to Type 1 rocks. FA700408 has a positive εNd value and does not plot close to any other group but may have a similar source to Southern Type 1 basalts with more crustal contamination. To illustrate mixing between Shaler Supergroup sedimentary rocks and Type 2 compositions in our study (represented by sample 93JP-71JB), mixing paths were calculated between the two endmembers with "x" marks denoting 10% intervals. Crustal contamination likely includes a small contribution from granitic basement rocks (not plotted) in addition to some sedimentary input (27).  (71). Gray data ellipses are excluded from the weighted mean dates. All analyses shown for samples S8 and FA700408 are from baddeleyite grains; all other samples only had zircon analyses (Figs. 3,4). For S8, the two analyses with gray data ellipses were excluded from the weighted mean but included in the upper-intercept age calculation. All uncertainties are 2σ. Uncertainties for weighted mean dates in this study are reported as ± X/Y/Z, where X represents internal error only, Y includes tracer calibration uncertainties, and Z includes both tracer calibration and decay constant uncertainties for comparisons with different isotopic chronometers. The reported uncertainties for the calculated upper-intercept dates for FA700408 and S8 represent the 95% confidence interval. 14RAT-513A is a coarse hornblende diorite that is mostly sericite-altered plagioclase and hornblende with chlorite alteration. B) S8 is a gabbro. It has intergranular textures with clinopyroxene growing between plagioclase crystals. C) 93JP-71JB is granitic with granophyric texture enclosing plagioclase crystals and secondary epidote. Sample 71M (not shown) is from lower in the same sill and is more gabbroic in composition, with intergranular textures of plagioclase and pyroxenes. D) 17RAT-R35B1, granite, also has granophyric textures of Kfeldspar and quartz around plagioclase crystals. E) 93JP-93K, hornblende quartz diorite with granophyric textures. F) 93JP-93L is also dioritic but more differentiated than 93JP-93K; the sample shares a similar mineralogy to 93JP-93K but has a greater proportion of feldspars and quartz and has extensive granophyric textures. G) FA700408 is a gabbro composed dominantly of plagioclase feldspar and clinopyroxene with minor granophyric K-feldspar and quartz. The clinopyroxene shows alteration to hornblende. H) F1966 is a hornblende diorite with granophyric texture.

Fig. S5. U-Th/He date vs. eU and t-T model paths.
Inverse thermal history modeling t-T path results for A) Model A: the Neoproterozoic Exhumation Hypothesis and B) Model B: the Neoproterozoic Burial/Exhumation Hypothesis. Color bar for acceptable fit paths is the same in all models. Black paths are good fit paths. White paths are representative good fit paths used to make predictions shown in C. Model constraints described in text and listed in Table S7 with all relevant model parameters. C) AHe (yellow) and ZHe (green) dates plotted against effective uranium concentration (eU). For zircon, grains with higher eU have higher levels of radiation damage and are sensitive to partial He loss at lower temperatures than lower eU grains. The two highest eU ZHe dates were excluded from modeling (see text). White symbols show forward model predictions of a representative path from each inverse model shown as white paths in A and B.      Table S5. U-Pb data table for all samples. Analyses in red were not included in calculations for mean sample dates. Analyses that were not run to completion due to high Pb blanks and/or low amounts of radiogenic Pb have been removed. The value 137.818 ± 0.045 (2σ) was used for the 238 U/ 235 U ratio in natural zircon (75). (a) z1, z2, etc. and b1, b2, etc. are labels for single zircon or baddeleyite grains/fragments, respectively. Zircon grains were annealed and chemically abraded after (16); (b) Model Th/U ratio iteratively calculated from the radiogenic 208 Pb/ 206 Pb ratio and 206 Pb/ 238 U age; (c) Pb* and Pbc represent radiogenic and common Pb, respectively; mol % 206 Pb* with respect to radiogenic, blank and initial common Pb; (d) Measured ratio corrected for spike and fractionation only. Fractionation estimated at 0.18 ± 0.03 %/a.m.u. for Daly analyses, based on analysis of NBS-981 and NBS-982; (e) Corrected for fractionation, spike, and common Pb; up to 1 pg of common Pb was assumed to be procedural blank: 206 Pb/ 204 Pb = 18.042 ± 0.61%; 207 Pb/ 204 Pb = 15.537 ± 0.52%; 208 Pb/ 204 Pb = 37.686 ± 0.63% (all uncertainties 1σ). Excess over blank was assigned to initial common Pb, using the (80) two-stage Pb isotope evolution model at the nominal sample age; (f) Errors are 2σ, propagated using the algorithms of (79); (g) Calculations are based on the decay constants of (100   (82). (b) Sample and mineral being analyzed. a is apatite. z is zircon. (c) Length is measured parallel to the c-axis and includes pyramidal terminations. It is measured twice on two perpendicular sides. (d) Width 1 is measured perpendicular to the c-axis. Width 2 is measured perpendicular to both the c-axis and width 1. (e) Geometry is defined as described as in Figure 3 of (88). 1 is ellipsoid, 2 is cylinder, 3 is tetrahedral prism, and 4 is hexagonal prism. f is noted if the analyzed grain is a fragment, otherwise the analyzed grain is a whole crystal. The combined alpha-ejection correction for the crystal calculated from the parent isotope-specific FT corrections, the proportion of U and Th contributing to 4 He production, and assuming homogeneous parent isotope distributions using equation A4 in (101). The parent isotope-specific alpha ejection-corrections were computed assuming the reported grain geometry in this table and the equations and alpha-stopping distances in (88). (u) The corrected (U-Th)/He date is calculated iteratively using the absolute values of He, U, Th and Sm, the isotope specific FT corrections, and equation 34 in (88) assuming secular equilibrium. (v) Uncertainty on the corrected (U-Th)/He date is reported at 2σ and includes the propagated total analytical uncertainties (TAU) on the U, Th, Sm and He measurements. Uncertainty propagation done using HeCalc (103). (x) Durango apatite fragments ran in conjunction with these analyses yield an unweighted mean and 2σ standard error of 30.8 ± 0.8 Ma (n=7). Fish Canyon Tuff zircon crystals ran in conjunction with these analyses yield an unweighted mean and 2σ standard error of 28.4 ± 1.0 Ma (n=7).